18  CJS [t,t] - growth model

Simple CJS model integrated with a growth in weight model to get phi, p, and growth estimates to compare with Neural Network CMR/growth model

In all the models below, 1 = not observed and 2 = observed in the input encounter data.
Encounter data are available here in the eh.csv file. Weight data are in weight.csv

18.1 modelNum 1: g(i,t) model

Growth rate (gr) model structure:

gr(t,i) <- grInt(t,i)

Model code is in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_makeFile.R
Not using this: Targets are loaded in R/modelCMR_tt_growth_NN_OB.R and saved as modelCMR_tt_growth_NN_OB_target

Model runs:

18.2 Retrieve model results

Code

library(targets)
# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html
# 
modelNum = 1 # growth only
dummy=2

#load('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_mostRecent.RData')
load(paste0('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_model', modelNum, '_mostRecent.RData'))
out_NN_OBmod1 <- d
rm(d)

#Input data
eh <- tar_read(eh_OB_2002_2014_target)

18.2.1 Model code

In the model code, a value of 1 for z or in gamma or omega signifies the individual is alive and a value of 2 signifies the individual is dead. y[i,j] is the observed data (encounter history file).

Code

out_NN_OBmod1$modelCode

{
    for (i in 1:N) {
        weight[i, first[i]] ~ dnorm(meanWeight_AIS[first[i]], 
            sd = 2)
        weightDATA[i, first[i]] ~ dnorm(weight[i, first[i]], 
            sd = weightSD)
        for (t in (first[i]):(last[i] - 1)) {
            weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i, 
                t]/90
            weightDATA[i, t + 1] ~ dnorm(weight[i, t + 1], sd = weightSD)
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
        }
    }
    for (t in 1:(T - 1)) {
        grIntT[t] ~ dnorm(0, sd = 1000)
    }
}

18.2.2 Model statistics

Run statistics

nIter	nBurnin	nChains	thinRate
5000	2500	2	5

[1] “Run time = 37.263 mins”

18.2.3 Plot traces

Code

  priors <- rnorm(out_NN_OBmod1$runData$nIter * out_NN_OBmod1$runData$nChains, 0, 1000)
  MCMCtrace(object = out_NN_OBmod1$mcmc,
            #ISB = FALSE,
            #exact = TRUE, 
            params = c("grIntT"),
            pdf = FALSE, 
            priors = priors)

18.2.4 Summary data

Code

  MCMCplot(object = out_NN_OBmod1$mcmc, params = c("grIntT"))

Summary values

Code

#|
  (summaryMod1_tt_growth <- MCMCsummary(object = out_NN_OBmod1$mcmc, params = c("grIntT"), round = 3) %>%
    rownames_to_column(var = "var"))

          var   mean    sd   2.5%    50%  97.5% Rhat n.eff
1   grIntT[1]  1.171 0.026  1.119  1.172  1.218 1.00   158
2   grIntT[2]  1.124 0.023  1.079  1.124  1.168 1.01    95
3   grIntT[3]  9.274 0.029  9.220  9.273  9.336 1.00    93
4   grIntT[4]  3.325 0.023  3.278  3.325  3.369 1.01   150
5   grIntT[5]  0.592 0.031  0.533  0.590  0.657 1.07   132
6   grIntT[6]  1.549 0.032  1.485  1.551  1.610 1.00   129
7   grIntT[7] 17.131 0.045 17.050 17.130 17.216 1.02    71
8   grIntT[8]  7.416 0.056  7.313  7.417  7.531 1.04    96
9   grIntT[9] -1.551 0.094 -1.717 -1.558 -1.356 1.00    35
10 grIntT[10]  4.074 0.095  3.875  4.078  4.254 1.03    32
11 grIntT[11] 22.904 0.115 22.682 22.903 23.116 1.34    34

Code

  # To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs
  source('./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R')
  
  #summaryMod1_tt_growth_z0 <- MCMCSummaryRMNA(object = out_NN_OBmod1$mcmc, params = c("weight"), round = 3) %>%
  summaryMod1_tt_growth_zGR <- MCMCSummaryRMNA(object = out_NN_OBmod1$mcmc, params = c("gr"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  summaryMod1_tt_growth_zW <- MCMCSummaryRMNA(object = out_NN_OBmod1$mcmc, params = c("weight"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  summaryMod1_tt_growth_z0 <- summaryMod1_tt_growth_zW |> 
    left_join(summaryMod1_tt_growth_zGR, by = c("ind", "t"), suffix = c("_weight", "_gr"))

Add other variables to summary values

Code

  ehLong <- 
    eh$eh |>
    as.data.frame() |> 
    rownames_to_column("ind0") |> 
    pivot_longer(starts_with("ais")) |> 
    mutate(
      t = as.numeric(str_match(name, "([0-9]+)")[,1]),
      enc = value,
      ind = as.numeric(ind0)
    ) |> 
    dplyr::select(-c(name, value, ind0))

  firstLong <- eh$first |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) |> 
    rename(first = value)
  
  lastLong <- eh$last |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) |> 
    rename(last = value)
    
  cohortLong <- eh$cohort |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) 
  
  weightLong <- 
    eh$weight |>
    as.data.frame() |> 
    rownames_to_column("ind0") |> 
    pivot_longer(starts_with("ais")) |> 
    mutate(
      t = as.numeric(str_match(name, "([0-9]+)")[,1]),
      weightDATA = value,
      ind = as.numeric(ind0)
    ) |> 
    dplyr::select(-c(name, value, ind0))
  
  summaryMod1_tt_growth_z <- summaryMod1_tt_growth_z0 |> 
    left_join(ehLong) |> 
    mutate(
      # meanM1 = mean - 1,
      # pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),
      enc8 = ifelse(enc == 1, 0.8, 0)
    ) |> 
    left_join(firstLong) |> 
    left_join(lastLong) |> 
    left_join(cohortLong) |> 
    left_join(weightLong) |> 
    mutate(
      resid = weightDATA - mean_weight
    ) 

Code

  resids <- summaryMod1_tt_growth_z |> 
    group_by(ind) |> 
    summarize(
      meanResid = mean(abs(resid), na.rm = TRUE),
      n = n()
    ) |> 
    ungroup()

  summaryMod1_tt_growth_zR <- 
    summaryMod1_tt_growth_z |> 
    left_join(resids)

Code

#  write.csv(summaryMod1_tt_growth_zR, './data/outForObservable/summaryMod1_tt_growth_zR.csv')

Code

ojs_define(numTags_OJS = dim(eh$tags)[1])
ojs_define(summaryMod1_tt_growth_zR_OJS = (summaryMod1_tt_growth_zR))

18.2.5 Plot predicted and observed masses for selected individuals

Select one or more individuals:

Code

viewof selectInd = Inputs.select([...Array(numTags_OJS).keys()], {
  label: "Which fish?",
  value: 1,
  step: 1,
  multiple: true
})

summaryMod1_tt_growth_zR_OJS_T = transpose(summaryMod1_tt_growth_zR_OJS)
summaryMod1_tt_growth_zR_OJS_T_selected = summaryMod1_tt_growth_zR_OJS_T.filter((d) =>
  selectInd.includes(d.ind)
)

Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort. Green dots are estimated growth rates.

Code

Plot.plot({
  width: width,
  height: 350,
  inset: 10,
  color: {
    scheme: "greys"
  },
  x: { label: "Age/season combination" },
  y: { label: "Body mass (mg)" },
  marks: [
    Plot.frame(),
    Plot.dot(summaryMod1_tt_growth_zR_OJS_T_selected, {
      x: "t",
      y: "mean_weight"
    }),
    Plot.line(summaryMod1_tt_growth_zR_OJS_T_selected, {
      x: "t",
      y: "mean_weight"
    }),
    Plot.dot(summaryMod1_tt_growth_zR_OJS_T_selected, {
      x: "t",
      y: "mean_gr",
      fill: "green"
    }),
    Plot.dot(summaryMod1_tt_growth_zR_OJS_T_selected, {
      x: "t",
      y: "weightDATA",
      fill: "orange"
    }),
    Plot.ruleX(summaryMod1_tt_growth_zR_OJS_T_selected, {
      x: "first",
      y: 3,
      stroke: "green"
    }),
    Plot.ruleX(summaryMod1_tt_growth_zR_OJS_T_selected, {
      x: "last",
      y: 3,
      stroke: "red"
    }),
    Plot.text(summaryMod1_tt_growth_zR_OJS_T_selected,
      Plot.selectFirst({
        x: 1,
        y: 60,
        frameAnchor: "top-left",
        text: (d) => "Cohort = " + d.cohort
      })
    ),
    Plot.text(summaryMod1_tt_growth_zR_OJS_T_selected,
      Plot.selectFirst({
        x: 1,
        y: 56,
        frameAnchor: "top-left",
        text: (d) => "Residual = " + d.meanResid
      })
    )
  ],
  facet: {
    data: summaryMod1_tt_growth_zR_OJS_T_selected,
    x: "ind"
  }
})

18.3 modelNum 2: phi(i,t) + g(i,t), p(i,t) model

Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)

Probability of survival (phi) model structure:

phi(t,i) <- phiInt(t,i)

Probability of capture (p) model structure:

p(t,i) <- pInt(t,i)

Growth rate (gr) model structure:

gr(t,i) <- grInt(t,i)

Survival/growth rate interaction:

Additive

Model code is in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_makeFile.R
Not using this: Targets are loaded in R/modelCMR_tt_growth_NN_OB.R and saved as modelCMR_tt_growth_NN_OB_target

Model runs:

18.3.1 Retrieve model results

Code

library(targets)
# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html
# 
modelNum = 2 # phi + growth

#load('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_mostRecent.RData')
load(paste0('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_model', modelNum, '_mostRecent.RData'))
out_NN_OBmod2 <- d
rm(d)

#Input data
eh <- tar_read(eh_OB_2002_2014_target)

18.3.2 Model code

In the model code, a value of 1 for z or in gamma or omega signifies the individual is alive and a value of 2 signifies the individual is dead. y[i,j] is the observed data (encounter history file).

Code

out_NN_OBmod2$modelCode

{
    for (i in 1:N) {
        weight[i, first[i]] ~ dnorm(meanWeight_AIS[first[i]], 
            sd = 2)
        weightDATA[i, first[i]] ~ dnorm(weight[i, first[i]], 
            sd = weightSD)
        for (t in (first[i]):(last[i] - 1)) {
            weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i, 
                t]/90
            weightDATA[i, t + 1] ~ dnorm(weight[i, t + 1], sd = weightSD)
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
        }
    }
    for (t in 1:(T - 1)) {
        grIntT[t] ~ dnorm(0, sd = 1000)
    }
    delta[1] <- 1
    delta[2] <- 0
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gamma[1, 1, t, i] <- phi[i, t]
            gamma[1, 2, t, i] <- 1 - phi[i, t]
            gamma[2, 1, t, i] <- 0
            gamma[2, 2, t, i] <- 1
        }
    }
    for (t in 1:(T - 1)) {
        p[t] ~ dunif(0, 1)
        omega[1, 1, t] <- 1 - p[t]
        omega[1, 2, t] <- p[t]
        omega[2, 1, t] <- 1
        omega[2, 2, t] <- 0
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            logit(phi[i, t]) <- phiInt[i, t]
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
        }
    }
    for (t in 1:(T - 1)) {
        phiIntT[t] ~ T(dnorm(0, sd = 1.5), -3.5, 3.5)
        phiIntTOut[t] <- ilogit(phiIntT[t])
    }
    for (i in 1:N) {
        z[i, first[i]] ~ dcat(delta[1:2])
        for (j in first[i]:(last[i] - 1)) {
            z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
            yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
        }
    }
}

18.3.3 Model statistics

Run statistics

nIter	nBurnin	nChains	thinRate
7500	5000	2	5

[1] “Run time = 2.231 hours”

18.3.4 Plot traces

Code

  priors <- rnorm(out_NN_OBmod2$runData$nIter * out_NN_OBmod2$runData$nChains, 0, 1000)
  MCMCtrace(object = out_NN_OBmod2$mcmc,
            #ISB = FALSE,
            #exact = TRUE, 
            params = c("phiIntTOut", "grIntT"),
            pdf = FALSE, 
            priors = priors)

18.3.5 Summary data

Code

  MCMCplot(object = out_NN_OBmod2$mcmc, params = c("phiIntTOut"))

Code

  MCMCplot(object = out_NN_OBmod2$mcmc, params = c("p"))

Code

  MCMCplot(object = out_NN_OBmod2$mcmc, params = c("grIntT"))

Summary values

Code

#|
(summaryMod2_tt_growth <- MCMCsummary(object = out_NN_OBmod2$mcmc, params = c("phiIntTOut", "p", "grIntT"), round = 3) %>%
    rownames_to_column(var = "var"))

              var   mean    sd   2.5%    50%  97.5% Rhat n.eff
1   phiIntTOut[1]  0.913 0.017  0.878  0.914  0.942 1.03    42
2   phiIntTOut[2]  0.903 0.017  0.867  0.905  0.930 1.12    55
3   phiIntTOut[3]  0.861 0.023  0.813  0.862  0.902 1.69    33
4   phiIntTOut[4]  0.729 0.034  0.655  0.731  0.799 2.26    39
5   phiIntTOut[5]  0.714 0.037  0.646  0.713  0.791 1.86    37
6   phiIntTOut[6]  0.791 0.036  0.703  0.797  0.845 1.10    32
7   phiIntTOut[7]  0.884 0.028  0.801  0.886  0.933 1.07    18
8   phiIntTOut[8]  0.258 0.040  0.177  0.260  0.330 1.03    25
9   phiIntTOut[9]  0.708 0.074  0.556  0.709  0.838 2.48    14
10 phiIntTOut[10]  0.675 0.086  0.484  0.686  0.850 2.59    13
11 phiIntTOut[11]  0.649 0.116  0.485  0.617  0.888 1.03    20
12           p[1]  0.595 0.021  0.557  0.595  0.637 1.01   282
13           p[2]  0.451 0.018  0.417  0.451  0.486 1.04   346
14           p[3]  0.706 0.021  0.668  0.706  0.747 1.30   126
15           p[4]  0.654 0.022  0.609  0.654  0.696 1.27   127
16           p[5]  0.633 0.025  0.584  0.633  0.680 1.05   220
17           p[6]  0.494 0.027  0.442  0.494  0.545 1.00   452
18           p[7]  0.730 0.035  0.663  0.729  0.798 1.01   113
19           p[8]  0.781 0.050  0.675  0.782  0.878 1.04   289
20           p[9]  0.634 0.062  0.513  0.637  0.749 1.31   439
21          p[10]  0.673 0.081  0.511  0.677  0.822 1.03   245
22          p[11]  0.802 0.125  0.552  0.815  0.987 1.07   102
23      grIntT[1]  1.172 0.028  1.117  1.171  1.229 1.06   119
24      grIntT[2]  1.120 0.023  1.076  1.119  1.168 1.01   106
25      grIntT[3]  9.281 0.029  9.222  9.281  9.337 1.20    92
26      grIntT[4]  3.324 0.022  3.281  3.324  3.366 1.00   281
27      grIntT[5]  0.583 0.030  0.525  0.583  0.638 1.00   117
28      grIntT[6]  1.558 0.032  1.495  1.557  1.620 1.09   139
29      grIntT[7] 17.114 0.041 17.035 17.113 17.200 1.10   102
30      grIntT[8]  7.416 0.055  7.306  7.417  7.528 1.02    81
31      grIntT[9] -1.559 0.083 -1.719 -1.559 -1.402 1.00    45
32     grIntT[10]  4.095 0.091  3.917  4.098  4.284 1.10    43
33     grIntT[11] 22.873 0.116 22.647 22.863 23.117 1.45    30

Code

  # To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs
  source('./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R')
  
  #summaryMod2_tt_growth_z0 <- MCMCSummaryRMNA(object = out_NN_OBmod2$mcmc, params = c("weight"), round = 3) %>%
  summaryMod2_tt_w <- MCMCSummaryRMNA(object = out_NN_OBmod2$mcmc, params = c("weight"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  summaryMod2_tt_z <- MCMCSummaryRMNA(object = out_NN_OBmod2$mcmc, params = c("z"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  summaryMod2_tt_z0 <- summaryMod2_tt_z |> 
    left_join(summaryMod2_tt_w, by = c("ind", "t"), suffix = c("_z", "_weight"))

Add other variables to summary values

Code

  ehLong <- 
    eh$eh |>
    as.data.frame() |> 
    rownames_to_column("ind0") |> 
    pivot_longer(starts_with("ais")) |> 
    mutate(
      t = as.numeric(str_match(name, "([0-9]+)")[,1]),
      enc = value,
      ind = as.numeric(ind0)
    ) |> 
    dplyr::select(-c(name, value, ind0))
  
  firstLong <- eh$first |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) |> 
    rename(first = value)
  
  lastLong <- eh$last |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) |> 
    rename(last = value)
    
  cohortLong <- eh$cohort |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) 
  
  weightLong <- 
    eh$weight |>
    as.data.frame() |> 
    rownames_to_column("ind0") |> 
    pivot_longer(starts_with("ais")) |> 
    mutate(
      t = as.numeric(str_match(name, "([0-9]+)")[,1]),
      weightDATA = value,
      ind = as.numeric(ind0)
    ) |> 
    dplyr::select(-c(name, value, ind0))
  
  summaryMod2_tt_z <- summaryMod2_tt_z0 |> 
    left_join(ehLong) |> 
    mutate(
      meanM1 = mean_z - 1,
      pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),
      enc8 = ifelse(enc == 1, 0.8, 0)
    ) |> 
    left_join(firstLong) |> 
    left_join(lastLong) |> 
    left_join(cohortLong) |> 
    left_join(weightLong) |> 
    mutate(
      resid = weightDATA - mean_weight
    ) 
    
  # summaryMod2_tt_z <- summaryMod2_tt_z0 |> 
  #   left_join(ehLong) |> 
  #   mutate(
  #     # meanM1 = mean - 1,
  #     # pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),
  #     enc8 = ifelse(enc == 1, 0.8, 0)
  #   ) |> 
  #   left_join(firstLong) |> 
  #   left_join(lastLong) |> 
  #   left_join(cohortLong) 

Code

  resids_mod2 <- summaryMod2_tt_z |> 
    group_by(ind) |> 
    summarize(
      meanResid = mean(abs(resid), na.rm = TRUE),
      n = n()
    ) |> 
    ungroup()

  summaryMod2_tt_zR <- 
    summaryMod2_tt_z |> 
    left_join(resids_mod2)

Code

ojs_define(numTags_OJS_mod2 = dim(eh$tags)[1])
ojs_define(summaryMod2_tt_zR_OJS = (summaryMod2_tt_zR))
#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod2))

18.3.6 Plot predicted and observed masses for selected individuals

Select one or more individuals:

Code

viewof selectInd_mod2 = Inputs.select([...Array(numTags_OJS).keys()], {
  label: "Which fish?",
  value: 1,
  step: 1,
  multiple: true
})

summaryMod2_tt_zR_OJS_T = transpose(summaryMod2_tt_zR_OJS)
summaryMod2_tt_zR_OJS_T_selected = summaryMod2_tt_zR_OJS_T.filter((d) =>
  selectInd_mod2.includes(d.ind)
)

Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Code

Plot.plot({
  width: width,
  height: 350,
  inset: 10,
  color: {
    scheme: "greys"
  },
  x: { label: "Age/season combination" },
  y: { label: "Body mass (mg)" },
  marks: [
    Plot.frame(),
    Plot.dot(summaryMod2_tt_zR_OJS_T_selected, {
      x: "t",
      y: "mean_weight"
    }),
    Plot.line(summaryMod2_tt_zR_OJS_T_selected, {
      x: "t",
      y: "mean_weight"
    }),
    <!-- Plot.dot(summaryMod2_tt_zR_OJS_T_selected, { -->
    <!--   x: "t", -->
    <!--   y: "mean_gr", -->
    <!--   fill: "green" -->
    <!-- }), -->
    Plot.dot(summaryMod2_tt_zR_OJS_T_selected, {
      x: "t",
      y: "weightDATA",
      fill: "orange"
    }),
    Plot.ruleX(summaryMod2_tt_zR_OJS_T_selected, {
      x: "first",
      y: 3,
      "stroke": "green"
    }),
    Plot.ruleX(summaryMod2_tt_zR_OJS_T_selected, {
      x: "last",
      y: 3,
      "stroke": "red"
    }),
    Plot.text(summaryMod2_tt_zR_OJS_T_selected,
      Plot.selectFirst({
        x: 1,
        y: 60,
        frameAnchor: "top-left",
        text: (d) => "Cohort = " + d.cohort
      })
    ),
    Plot.text(summaryMod2_tt_zR_OJS_T_selected,
      Plot.selectFirst({
        x: 1,
        y: 56,
        frameAnchor: "top-left",
        text: (d) => "Residual = " + d.meanResid
      })
    )
  ],
  facet: {
    data: summaryMod2_tt_zR_OJS_T_selected,
    x: "ind"
  }
})

18.3.7 Plot survival

Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Code

Plot.plot({
  width: width,
  height: 350,
  inset: 10,
  color: {
    scheme: "greys"
  },
  x: { label: "Age/season combination" },
  y: { label: "Probability of survival" },
  marks: [
    Plot.frame(),
    Plot.dot(summaryMod2_tt_zR_OJS_T_selected, {
      x: "t",
      y: "pSurv"
    }),
    Plot.dot(summaryMod2_tt_zR_OJS_T_selected, {
      x: "t",
      y: "enc8",
      fill: "orange"
    }),
    Plot.line(summaryMod2_tt_zR_OJS_T_selected, {
      x: "t",
      y: "pSurv"
    }),
    Plot.ruleX(summaryMod2_tt_zR_OJS_T_selected, {
      x: "first",
      y: 1,
      "stroke": "green"
    }),
    Plot.ruleX(summaryMod2_tt_zR_OJS_T_selected, {
      x: "last",
      y: 1,
      "stroke": "red"
    })
    ,
    Plot.text(summaryMod2_tt_zR_OJS_T_selected, 
      Plot.selectFirst({
        x: 11,
        y: 1,
        text: d => d.cohort
      })
    )
  ],
  facet: {
    data: summaryMod2_tt_zR_OJS_T_selected,
    x: "ind"
  }
})

18.4 modelNum 3: phi(i,t) * g(i,t), p(i,t) model

Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)

Probability of survival (phi) model structure:

phi(t,i) <- phiInt(t,i) + weight(t,i) + weight(t,i)^2

Probability of capture (p) model structure:

p(t,i) <- pInt(t,i)

Growth rate (gr) model structure:

gr(t,i) <- grInt(t,i)

Survival/growth rate interaction:

Multiplicative

Model code is in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_makeFile.R
Not using this: Targets are loaded in R/modelCMR_tt_growth_NN_OB.R and saved as modelCMR_tt_growth_NN_OB_target

Model runs:

18.4.1 Retrieve model results

Code

library(targets)
# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html
# 
modelNum = 3 # phi + growth

#load('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_mostRecent.RData')
load(paste0('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_model', modelNum, '_mostRecent.RData'))
out_NN_OBmod3 <- d
rm(d)

#Input data
eh <- tar_read(eh_OB_2002_2014_target)

18.4.2 Model code

In the model code, a value of 1 for z or in gamma or omega signifies the individual is alive and a value of 2 signifies the individual is dead. y[i,j] is the observed data (encounter history file).

Code

out_NN_OBmod3$modelCode

{
    for (i in 1:N) {
        weight[i, first[i]] ~ dnorm(meanWeight_AISstd[first[i]], 
            sd = 2)
        weightDATAstd[i, first[i]] ~ dnorm(weight[i, first[i]], 
            sd = weightSD)
        for (t in (first[i]):(last[i] - 1)) {
            weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i, 
                t]/90
            weightDATAstd[i, t + 1] ~ dnorm(weight[i, t + 1], 
                sd = weightSD)
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
        }
    }
    for (t in 1:(T - 1)) {
        grIntT[t] ~ dnorm(0, sd = 1000)
    }
    delta[1] <- 1
    delta[2] <- 0
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gamma[1, 1, t, i] <- phi[i, t]
            gamma[1, 2, t, i] <- 1 - phi[i, t]
            gamma[2, 1, t, i] <- 0
            gamma[2, 2, t, i] <- 1
        }
    }
    for (t in 1:(T - 1)) {
        p[t] ~ dunif(0, 1)
        omega[1, 1, t] <- 1 - p[t]
        omega[1, 2, t] <- p[t]
        omega[2, 1, t] <- 1
        omega[2, 2, t] <- 0
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            logit(phi[i, t]) <- phiInt[i, t] + phiBeta[1, i, 
                t] * weight[i, t] + phiBeta[2, i, t] * weight[i, 
                t] * weight[i, t]
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
            for (b in 1:2) {
                phiBeta[b, i, t] ~ dnorm(phiBetaT[b, t], sd = 1/0.667)
            }
        }
    }
    for (t in 1:(T - 1)) {
        phiIntT[t] ~ T(dnorm(0, sd = 0.667), -3.5, 3.5)
        phiIntTOut[t] <- ilogit(phiIntT[t])
        phiBetaT[1, t] ~ dnorm(0, sd = 1/0.667)
        phiBetaT[2, t] ~ dnorm(0, sd = 1/0.667)
    }
    for (i in 1:N) {
        z[i, first[i]] ~ dcat(delta[1:2])
        for (j in first[i]:(last[i] - 1)) {
            z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
            yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
        }
    }
}

18.4.3 Model statistics

Run statistics

nIter	nBurnin	nChains	thinRate
15000	10000	2	5

[1] “Run time = 6.035 hours”

18.4.4 Plot traces

Code

  priors <- rnorm(out_NN_OBmod3$runData$nIter * out_NN_OBmod3$runData$nChains, 0, 1000)
  MCMCtrace(object = out_NN_OBmod3$mcmc,
            #ISB = FALSE,
            #exact = TRUE, 
            params = c("phiIntTOut", "grIntT"),
            pdf = FALSE, 
            priors = priors)

18.4.5 Summary data

Code

  MCMCplot(object = out_NN_OBmod3$mcmc, params = c("phiIntTOut"))

Code

  MCMCplot(object = out_NN_OBmod3$mcmc, params = c("phiBetaT"))

Code

  MCMCplot(object = out_NN_OBmod3$mcmc, params = c("p"))

Code

  MCMCplot(object = out_NN_OBmod3$mcmc, params = c("grIntT"))

Summary values

Code

#|
(summaryMod3_tt_growth <- MCMCsummary(object = out_NN_OBmod3$mcmc, params = c("phiIntTOut", "p", "grIntT", "phiBetaT"), round = 3) %>%
    rownames_to_column(var = "var"))

               var   mean    sd   2.5%    50%  97.5% Rhat n.eff
1    phiIntTOut[1]  0.809 0.063  0.674  0.820  0.906 1.00    20
2    phiIntTOut[2]  0.882 0.047  0.788  0.887  0.951 4.70    30
3    phiIntTOut[3]  0.795 0.046  0.696  0.797  0.875 2.13    26
4    phiIntTOut[4]  0.748 0.031  0.680  0.749  0.802 1.40    82
5    phiIntTOut[5]  0.755 0.032  0.686  0.755  0.813 1.06    74
6    phiIntTOut[6]  0.809 0.038  0.735  0.809  0.876 1.22    33
7    phiIntTOut[7]  0.867 0.035  0.781  0.871  0.918 1.19    31
8    phiIntTOut[8]  0.470 0.092  0.293  0.463  0.635 2.15    22
9    phiIntTOut[9]  0.592 0.095  0.394  0.590  0.794 1.15    16
10  phiIntTOut[10]  0.596 0.112  0.338  0.599  0.813 1.58    16
11  phiIntTOut[11]  0.647 0.113  0.421  0.669  0.823 1.24    10
12            p[1]  0.590 0.021  0.549  0.590  0.630 1.04   177
13            p[2]  0.409 0.018  0.374  0.408  0.446 1.05   104
14            p[3]  0.704 0.019  0.665  0.704  0.741 1.01   306
15            p[4]  0.655 0.021  0.610  0.655  0.696 1.01   323
16            p[5]  0.635 0.024  0.587  0.636  0.684 1.00   311
17            p[6]  0.495 0.028  0.440  0.495  0.549 1.01   399
18            p[7]  0.746 0.033  0.682  0.745  0.808 1.06   318
19            p[8]  0.774 0.049  0.672  0.775  0.863 1.16   206
20            p[9]  0.596 0.067  0.459  0.597  0.725 1.24    79
21           p[10]  0.657 0.103  0.446  0.668  0.834 1.22    30
22           p[11]  0.812 0.133  0.531  0.832  0.993 1.05    24
23       grIntT[1]  0.122 0.027  0.064  0.122  0.171 1.01   292
24       grIntT[2]  0.128 0.024  0.081  0.128  0.178 1.17   175
25       grIntT[3]  0.899 0.029  0.842  0.899  0.956 1.35   166
26       grIntT[4]  0.327 0.024  0.282  0.325  0.377 1.12   226
27       grIntT[5]  0.075 0.033  0.011  0.075  0.140 1.03   140
28       grIntT[6]  0.170 0.033  0.101  0.171  0.234 1.04   170
29       grIntT[7]  1.642 0.044  1.550  1.643  1.724 1.02   117
30       grIntT[8]  0.722 0.053  0.616  0.719  0.826 1.02    80
31       grIntT[9] -0.164 0.091 -0.360 -0.165  0.014 1.08    28
32      grIntT[10]  0.435 0.104  0.224  0.443  0.655 1.08    25
33      grIntT[11]  2.150 0.119  1.924  2.152  2.385 1.13    41
34  phiBetaT[1, 1] -2.174 0.657 -3.506 -2.077 -1.092 1.13    12
35  phiBetaT[2, 1] -0.324 0.658 -1.536 -0.377  1.361 1.18    13
36  phiBetaT[1, 2] -1.027 0.970 -2.812 -1.071  0.732 5.66     7
37  phiBetaT[2, 2]  1.928 0.568  0.778  1.887  3.083 1.15    11
38  phiBetaT[1, 3] -3.166 0.528 -3.974 -3.276 -2.059 3.41    12
39  phiBetaT[2, 3] -2.785 0.549 -3.711 -2.847 -1.520 1.14    11
40  phiBetaT[1, 4] -0.430 0.429 -1.190 -0.465  0.477 1.80    15
41  phiBetaT[2, 4] -0.295 0.488 -1.244 -0.340  0.735 3.06    16
42  phiBetaT[1, 5] -0.824 0.337 -1.427 -0.849 -0.151 1.60    19
43  phiBetaT[2, 5]  0.176 0.349 -0.459  0.189  0.863 1.52    19
44  phiBetaT[1, 6] -1.562 0.442 -2.438 -1.608 -0.550 1.98    13
45  phiBetaT[2, 6]  1.016 0.417  0.263  1.035  1.728 3.74    26
46  phiBetaT[1, 7] -0.577 0.832 -2.419 -0.659  0.896 2.99     6
47  phiBetaT[2, 7]  0.856 0.583 -0.120  0.905  1.966 2.63    11
48  phiBetaT[1, 8] -0.664 0.670 -1.946 -0.718  0.485 5.57    10
49  phiBetaT[2, 8] -0.421 0.356 -1.007 -0.438  0.276 3.74    21
50  phiBetaT[1, 9]  1.117 0.890 -0.927  1.363  2.444 1.00     6
51  phiBetaT[2, 9]  0.277 0.547 -0.740  0.228  1.338 1.63    11
52 phiBetaT[1, 10]  0.475 0.519 -0.451  0.476  1.442 1.30    12
53 phiBetaT[2, 10]  0.052 0.542 -1.060  0.074  1.002 1.77     9
54 phiBetaT[1, 11]  0.583 0.730 -0.578  0.474  1.929 1.91     6
55 phiBetaT[2, 11] -0.083 0.480 -1.066 -0.032  0.828 1.44    13

Code

  # To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs
  source('./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R')
  
  #summaryMod3_tt_growth_z0 <- MCMCSummaryRMNA(object = out_NN_OBmod3$mcmc, params = c("weight"), round = 3) %>%
  summaryMod3_tt_w <- MCMCSummaryRMNA(object = out_NN_OBmod3$mcmc, params = c("weight"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  summaryMod3_tt_z <- MCMCSummaryRMNA(object = out_NN_OBmod3$mcmc, params = c("z"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  summaryMod3_tt_PhiBetaT <- MCMCSummaryRMNA(object = out_NN_OBmod3$mcmc, params = c("phiBetaT"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  summaryMod3_tt_z0 <- summaryMod3_tt_z |> 
    left_join(summaryMod3_tt_w, by = c("ind", "t"), suffix = c("_z", "_weight"))

Code

summaryMod3_tt_growth_phi <- MCMCsummary(object = out_NN_OBmod3$mcmc, params = c("phiIntTOut", "grIntT"), round = 3) |> 
    rownames_to_column(var = "var") |> 
    mutate(
      t = as.numeric(str_match(var, "[0-9]+"))
    ) |> 
    separate_wider_delim(var, "[", names = c("var1", "var2")) |> 
  dplyr::select(mean,var1,t) |> 
  pivot_wider(names_from = var1, names_prefix = "beta", values_from = mean)

summaryMod3_tt_growth_phiBetas = MCMCsummary(object = out_NN_OBmod3$mcmc, params = c("phiBetaT"), round = 3) |> 
        rownames_to_column(var = "var") |> 
        mutate(
          beta = str_match(var, "([0-9]+), ([0-9]+)")[,2],
          t = as.numeric(str_match(var, "([0-9]+), ([0-9]+)")[,3])
        ) |> 
  dplyr::select(mean,beta,t) |> 
  pivot_wider(names_from = beta, names_prefix = "beta", values_from = mean)



summaryMod3_tt_growth_phiParams <-  
  bind_cols(summaryMod3_tt_growth_phi, summaryMod3_tt_growth_phiBetas |> dplyr::select(-t)) |> 
  mutate(
    phiInt = logit(betaphiIntTOut)
  )

weights = expand.grid(weights = seq(-2,2,0.25), t = 1:11)
#add predictions...

Add other variables to summary values

Code

  ehLong <- 
    eh$eh |>
    as.data.frame() |> 
    rownames_to_column("ind0") |> 
    pivot_longer(starts_with("ais")) |> 
    mutate(
      t = as.numeric(str_match(name, "([0-9]+)")[,1]),
      enc = value,
      ind = as.numeric(ind0)
    ) |> 
    dplyr::select(-c(name, value, ind0))
  
  firstLong <- eh$first |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) |> 
    rename(first = value)
  
  lastLong <- eh$last |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) |> 
    rename(last = value)
    
  cohortLong <- eh$cohort |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) 
  
  weightLong <- 
    eh$weight |>
    as.data.frame() |> 
    rownames_to_column("ind0") |> 
    pivot_longer(starts_with("ais")) |> 
    mutate(
      t = as.numeric(str_match(name, "([0-9]+)")[,1]),
      weightDATA = value,
      ind = as.numeric(ind0)
    ) |> 
    dplyr::select(-c(name, value, ind0))
  
  meanWeight <- mean(eh$weight, na.rm = TRUE)
  sdWeight <- sd(eh$weight, na.rm = TRUE)

  summaryMod3_tt_z <- summaryMod3_tt_z0 |> 
    left_join(ehLong) |> 
    mutate(
      meanM1 = mean_z - 1,
      pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),
      enc8 = ifelse(enc == 1, 0.8, 0)
    ) |> 
    left_join(firstLong) |> 
    left_join(lastLong) |> 
    left_join(cohortLong) |> 
    left_join(weightLong) |> 
    mutate(
      weightDATA_std = (weightDATA - meanWeight) / sdWeight, 
      resid = weightDATA_std - mean_weight
    ) 
    
  # summaryMod3_tt_z <- summaryMod3_tt_z0 |> 
  #   left_join(ehLong) |> 
  #   mutate(
  #     # meanM1 = mean - 1,
  #     # pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),
  #     enc8 = ifelse(enc == 1, 0.8, 0)
  #   ) |> 
  #   left_join(firstLong) |> 
  #   left_join(lastLong) |> 
  #   left_join(cohortLong) 

Code

  resids_mod3 <- summaryMod3_tt_z |> 
    group_by(ind) |> 
    summarize(
      meanResid = mean(abs(resid), na.rm = TRUE),
      n = n()
    ) |> 
    ungroup()

  summaryMod3_tt_zR <- 
    summaryMod3_tt_z |> 
    left_join(resids_mod3)

Code

ojs_define(numTags_OJS_mod3 = dim(eh$tags)[1])
ojs_define(summaryMod3_tt_zR_OJS = (summaryMod3_tt_zR))
#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod3))

18.4.6 Plot predicted and observed masses for selected individuals

Select one or more individuals:

Code

viewof selectInd_mod3 = Inputs.select([...Array(numTags_OJS).keys()], {
  label: "Which fish?",
  value: 1,
  step: 1,
  multiple: true
})

summaryMod3_tt_zR_OJS_T = transpose(summaryMod3_tt_zR_OJS)
summaryMod3_tt_zR_OJS_T_selected = summaryMod3_tt_zR_OJS_T.filter((d) =>
  selectInd_mod3.includes(d.ind)
)

Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Code

Plot.plot({
  width: width,
  height: 350,
  inset: 10,
  color: {
    scheme: "greys"
  },
  x: { label: "Age/season combination" },
  y: { label: "Standardized body mass (mg)" },
  marks: [
    Plot.frame(),
    Plot.dot(summaryMod3_tt_zR_OJS_T_selected, {
      x: "t",
      y: "mean_weight"
    }),
    Plot.line(summaryMod3_tt_zR_OJS_T_selected, {
      x: "t",
      y: "mean_weight"
    }),
    <!-- Plot.dot(summaryMod3_tt_zR_OJS_T_selected, { -->
    <!--   x: "t", -->
    <!--   y: "mean_gr", -->
    <!--   fill: "green" -->
    <!-- }), -->
    Plot.dot(summaryMod3_tt_zR_OJS_T_selected, {
      x: "t",
      y: "weightDATA_std",
      fill: "orange"
    }),
    Plot.ruleX(summaryMod3_tt_zR_OJS_T_selected, {
      x: "first",
      y: 3,
      "stroke": "green"
    }),
    Plot.ruleX(summaryMod3_tt_zR_OJS_T_selected, {
      x: "last",
      y: 3,
      "stroke": "red"
    }),
    Plot.text(summaryMod3_tt_zR_OJS_T_selected,
      Plot.selectFirst({
        x: 1,
        y: 4,
        frameAnchor: "top-left",
        text: (d) => "Cohort = " + d.cohort
      })
    ),
    Plot.text(summaryMod3_tt_zR_OJS_T_selected,
      Plot.selectFirst({
        x: 1,
        y: 3.5,
        frameAnchor: "top-left",
        text: (d) => "Residual = " + d.meanResid
      })
    )
  ],
  facet: {
    data: summaryMod3_tt_zR_OJS_T_selected,
    x: "ind"
  }
})

18.4.7 Plot survival

Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Code

Plot.plot({
  width: width,
  height: 350,
  inset: 10,
  color: {
    scheme: "greys"
  },
  x: { label: "Age/season combination" },
  y: { label: "Probability of survival" },
  marks: [
    Plot.frame(),
    Plot.dot(summaryMod3_tt_zR_OJS_T_selected, {
      x: "t",
      y: "pSurv"
    }),
    Plot.dot(summaryMod3_tt_zR_OJS_T_selected, {
      x: "t",
      y: "enc8",
      fill: "orange"
    }),
    Plot.line(summaryMod3_tt_zR_OJS_T_selected, {
      x: "t",
      y: "pSurv"
    }),
    Plot.ruleX(summaryMod3_tt_zR_OJS_T_selected, {
      x: "first",
      y: 1,
      "stroke": "green"
    }),
    Plot.ruleX(summaryMod3_tt_zR_OJS_T_selected, {
      x: "last",
      y: 1,
      "stroke": "red"
    })
    ,
    Plot.text(summaryMod3_tt_zR_OJS_T_selected, 
      Plot.selectFirst({
        x: 11,
        y: 1,
        text: d => d.cohort
      })
    )
  ],
  facet: {
    data: summaryMod3_tt_zR_OJS_T_selected,
    x: "ind"
  }
})

18.5 modelNum 4: phi(i,ys) * g(i,ys), p(i,t) model

Model with phi and p for each yoy(y)/season(s) combination and individual (i)

Probability of survival (phi) model structure:

phi(t,i) <- phiInt(ys,i) + weight(t,i) + weight(t,i)^2

Probability of capture (p) model structure:

p(t,i) <- pInt(t,i)

Growth rate (gr) model structure:

gr(t,i) <- grInt(ys,i)

Survival/growth rate interaction:

Multiplicative

Model code is in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_makeFile.R
Not using this: Targets are loaded in R/modelCMR_tt_growth_NN_OB.R and saved as modelCMR_tt_growth_NN_OB_target

Model runs:

18.5.1 Retrieve model results

Code

library(targets)
# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html
# 
modelNum = 4 # phi + growth

#load('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_mostRecent.RData')
load(paste0('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_model', modelNum, '_mostRecent.RData'))
out_NN_OBmod4 <- d
rm(d)

#Input data
eh <- tar_read(eh_OB_2002_2014_target)

18.5.2 Model code

In the model code, a value of 1 for z or in gamma or omega signifies the individual is alive and a value of 2 signifies the individual is dead. y[i,j] is the observed data (encounter history file).

Code

out_NN_OBmod4$modelCode

{
    for (i in 1:N) {
        weight[i, first[i]] ~ dnorm(meanWeight_AISstd[first[i]], 
            sd = 2)
        weightDATAstd[i, first[i]] ~ dnorm(weight[i, first[i]], 
            sd = weightSD)
        for (t in (first[i]):(last[i] - 1)) {
            weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i, 
                t]/90
            weightDATAstd[i, t + 1] ~ dnorm(weight[i, t + 1], 
                sd = weightSD)
            gr[i, t] <- grInt[isYOY_DATA[i, t], seasons[t]]
        }
    }
    for (yoy in 1:2) {
        for (s in 1:4) {
            grInt[yoy, s] ~ dnorm(0, sd = 100)
        }
    }
    delta[1] <- 1
    delta[2] <- 0
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gamma[1, 1, t, i] <- phi[i, t]
            gamma[1, 2, t, i] <- 1 - phi[i, t]
            gamma[2, 1, t, i] <- 0
            gamma[2, 2, t, i] <- 1
        }
    }
    for (t in 1:(T - 1)) {
        p[t] ~ dunif(0, 1)
        omega[1, 1, t] <- 1 - p[t]
        omega[1, 2, t] <- p[t]
        omega[2, 1, t] <- 1
        omega[2, 2, t] <- 0
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            logit(phi[i, t]) <- phiInt[isYOY_DATA[i, t], seasons[t]] + 
                phiBeta[1, isYOY_DATA[i, t], seasons[t]] * weight[i, 
                  t] + phiBeta[2, isYOY_DATA[i, t], seasons[t]] * 
                weight[i, t] * weight[i, t]
        }
    }
    for (yoy in 1:2) {
        for (s in 1:4) {
            phiInt[yoy, s] ~ T(dnorm(0, sd = 1/0.667), -3.5, 
                3.5)
            for (b in 1:2) {
                phiBeta[b, yoy, s] ~ dnorm(0, sd = 1/0.667)
            }
        }
    }
    for (i in 1:N) {
        z[i, first[i]] ~ dcat(delta[1:2])
        for (j in first[i]:(last[i] - 1)) {
            z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
            yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
        }
    }
}

18.5.3 Model statistics

Run statistics

nIter	nBurnin	nChains	thinRate
5000	2500	2	5

[1] “Run time = 1.59 hours”

18.5.4 Plot traces

Code

  priors <- rnorm(out_NN_OBmod4$runData$nIter * out_NN_OBmod4$runData$nChains, 0, 1000)
  MCMCtrace(object = out_NN_OBmod4$mcmc,
            #ISB = FALSE,
            #exact = TRUE, 
            params = c("phiInt", "grIntYS"),
            pdf = FALSE, 
            priors = priors)

18.5.5 Summary data

Code

  MCMCplot(object = out_NN_OBmod4$mcmc, params = c("phiInt"))

Code

  MCMCplot(object = out_NN_OBmod4$mcmc, params = c("phiBeta"))

Code

  MCMCplot(object = out_NN_OBmod4$mcmc, params = c("p"))

Code

  MCMCplot(object = out_NN_OBmod4$mcmc, params = c("grInt"))

Summary values

Code

#|
(summaryMod4_tt_growth <- MCMCsummary(object = out_NN_OBmod4$mcmc, params = c("phiInt", "p", "grInt", "phiBeta"), round = 3) %>%
    rownames_to_column(var = "var"))

                var   mean     sd     2.5%    50%   97.5% Rhat n.eff
1      phiInt[1, 1]  2.185  0.301    1.662  2.154   2.849 1.03    82
2      phiInt[2, 1]  3.070  0.303    2.366  3.123   3.479 1.04   213
3      phiInt[1, 2] -0.027  1.411   -2.502 -0.058   2.789 1.01  1000
4      phiInt[2, 2]  3.148  0.208    2.687  3.166   3.467 1.12    68
5      phiInt[1, 3]  3.002  0.327    2.290  3.031   3.459 1.00   151
6      phiInt[2, 3]  2.819  0.349    2.142  2.845   3.416 1.06    81
7      phiInt[1, 4]  2.923  0.331    2.279  2.950   3.459 1.03   123
8      phiInt[2, 4]  2.620  0.342    1.899  2.608   3.282 1.16    73
9              p[1]  0.508  0.018    0.473  0.507   0.541 1.01   595
10             p[2]  0.342  0.014    0.315  0.343   0.373 1.01   326
11             p[3]  0.465  0.014    0.439  0.465   0.493 1.01   579
12             p[4]  0.328  0.012    0.304  0.328   0.350 1.01   756
13             p[5]  0.225  0.010    0.205  0.225   0.245 1.02   621
14             p[6]  0.135  0.008    0.120  0.135   0.151 1.10   820
15             p[7]  0.164  0.009    0.146  0.164   0.183 1.02   762
16             p[8]  0.067  0.006    0.055  0.067   0.080 1.02   734
17             p[9]  0.036  0.005    0.027  0.036   0.045 1.02   840
18            p[10]  0.024  0.004    0.017  0.023   0.032 1.02   686
19            p[11]  0.018  0.003    0.012  0.018   0.025 1.01   780
20      grInt[1, 1]  0.839  0.015    0.811  0.838   0.870 1.00   136
21      grInt[2, 1]  1.802  0.022    1.757  1.802   1.845 1.02   201
22      grInt[1, 2]  1.965 98.894 -184.213  3.215 200.363 1.00  1055
23      grInt[2, 2]  0.391  0.009    0.374  0.390   0.409 1.02   187
24      grInt[1, 3]  0.097  0.013    0.069  0.098   0.123 1.00   180
25      grInt[2, 3] -0.043  0.015   -0.072 -0.043  -0.016 1.01   189
26      grInt[1, 4]  0.159  0.009    0.140  0.158   0.177 1.02   141
27      grInt[2, 4]  0.202  0.012    0.179  0.202   0.225 1.03   148
28 phiBeta[1, 1, 1]  1.349  0.712    0.042  1.383   2.745 1.14    49
29 phiBeta[2, 1, 1]  3.134  0.925    1.475  3.095   5.143 1.27    70
30 phiBeta[1, 2, 1]  0.194  0.874   -1.401  0.150   2.010 1.11   165
31 phiBeta[2, 2, 1]  2.021  0.940    0.465  1.956   4.005 1.00   160
32 phiBeta[1, 1, 2] -0.034  1.531   -3.003 -0.043   3.158 1.00  1084
33 phiBeta[2, 1, 2] -0.008  1.492   -2.802  0.025   2.938 1.00   927
34 phiBeta[1, 2, 2] -0.654  0.372   -1.314 -0.660   0.027 3.71    21
35 phiBeta[2, 2, 2]  0.153  0.126   -0.077  0.147   0.390 3.30    25
36 phiBeta[1, 1, 3]  0.548  0.929   -1.147  0.507   2.429 1.09    47
37 phiBeta[2, 1, 3]  0.909  1.073   -0.877  0.887   3.453 1.13    38
38 phiBeta[1, 2, 3] -0.410  0.865   -2.644 -0.303   1.005 1.07    43
39 phiBeta[2, 2, 3]  1.576  0.943    0.123  1.526   3.559 1.00    32
40 phiBeta[1, 1, 4]  1.127  0.986   -0.633  1.123   3.444 1.16    25
41 phiBeta[2, 1, 4]  1.857  0.917    0.302  1.864   3.821 1.14    65
42 phiBeta[1, 2, 4] -0.309  0.654   -1.639 -0.301   1.019 1.01   100
43 phiBeta[2, 2, 4]  1.797  0.937    0.449  1.695   3.866 1.02    56

Code

  # To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs
  source('./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R')
  
  #summaryMod4_tt_growth_z0 <- MCMCSummaryRMNA(object = out_NN_OBmod4$mcmc, params = c("weight"), round = 4) %>%
  summaryMod4_tt_w <- MCMCSummaryRMNA(object = out_NN_OBmod4$mcmc, params = c("weight"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  summaryMod4_tt_z <- MCMCSummaryRMNA(object = out_NN_OBmod4$mcmc, params = c("z"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  # summaryMod4_tt_PhiBetaT <- MCMCSummaryRMNA(object = out_NN_OBmod4$mcmc, params = c("phiBeta"), round = 3) %>%
  #   rownames_to_column(var = "var") |> 
  #   mutate(
  #     ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
  #     t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
  #     ind = as.numeric(ind0),
  #     t = as.numeric(t0)
  #   ) |> 
  #   dplyr::select(-c(t0, ind0))
  
  summaryMod4_tt_z0 <- summaryMod4_tt_z |> 
    left_join(summaryMod4_tt_w, by = c("ind", "t"), suffix = c("_z", "_weight"))

Add other variables to summary values

Code

  ehLong <- 
    eh$eh |>
    as.data.frame() |> 
    rownames_to_column("ind0") |> 
    pivot_longer(starts_with("ais")) |> 
    mutate(
      t = as.numeric(str_match(name, "([0-9]+)")[,1]),
      enc = value,
      ind = as.numeric(ind0)
    ) |> 
    dplyr::select(-c(name, value, ind0))
  
  firstLong <- eh$first |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) |> 
    rename(first = value)
  
  lastLong <- eh$last |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) |> 
    rename(last = value)
    
  cohortLong <- eh$cohort |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) 
  
  weightLong <- 
    eh$weight |>
    as.data.frame() |> 
    rownames_to_column("ind0") |> 
    pivot_longer(starts_with("ais")) |> 
    mutate(
      t = as.numeric(str_match(name, "([0-9]+)")[,1]),
      weightDATA = value,
      ind = as.numeric(ind0)
    ) |> 
    dplyr::select(-c(name, value, ind0))
  
  meanWeight <- mean(eh$weight, na.rm = TRUE)
  sdWeight <- sd(eh$weight, na.rm = TRUE)

  summaryMod4_tt_z <- summaryMod4_tt_z0 |> 
    left_join(ehLong) |> 
    mutate(
      meanM1 = mean_z - 1,
      pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),
      enc8 = ifelse(enc == 1, 0.8, 0)
    ) |> 
    left_join(firstLong) |> 
    left_join(lastLong) |> 
    left_join(cohortLong) |> 
    left_join(weightLong) |> 
    mutate(
      weightDATA_std = (weightDATA - meanWeight) / sdWeight, 
      resid = weightDATA_std - mean_weight
    ) 
    
  # summaryMod4_tt_z <- summaryMod4_tt_z0 |> 
  #   left_join(ehLong) |> 
  #   mutate(
  #     # meanM1 = mean - 1,
  #     # pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),
  #     enc8 = ifelse(enc == 1, 0.8, 0)
  #   ) |> 
  #   left_join(firstLong) |> 
  #   left_join(lastLong) |> 
  #   left_join(cohortLong) 

Code

  resids_mod4 <- summaryMod4_tt_z |> 
    group_by(ind) |> 
    summarize(
      meanResid = mean(abs(resid), na.rm = TRUE),
      n = n()
    ) |> 
    ungroup()

  summaryMod4_tt_zR <- 
    summaryMod4_tt_z |> 
    left_join(resids_mod4)

Code

ojs_define(numTags_OJS_mod4 = dim(eh$tags)[1])
ojs_define(summaryMod4_tt_zR_OJS = (summaryMod4_tt_zR))
#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod4))

18.5.6 Plot predicted and observed masses for selected individuals

Select one or more individuals:

Code

viewof selectInd_mod4 = Inputs.select([...Array(numTags_OJS).keys()], {
  label: "Which fish?",
  value: 1,
  step: 1,
  multiple: true
})

summaryMod4_tt_zR_OJS_T = transpose(summaryMod4_tt_zR_OJS)
summaryMod4_tt_zR_OJS_T_selected = summaryMod4_tt_zR_OJS_T.filter((d) =>
  selectInd_mod4.includes(d.ind)
)

Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Code

Plot.plot({
  width: width,
  height: 350,
  inset: 10,
  color: {
    scheme: "greys"
  },
  x: { label: "Age/season combination" },
  y: { label: "Standardized body mass (mg)" },
  marks: [
    Plot.frame(),
    Plot.dot(summaryMod4_tt_zR_OJS_T_selected, {
      x: "t",
      y: "mean_weight"
    }),
    Plot.line(summaryMod4_tt_zR_OJS_T_selected, {
      x: "t",
      y: "mean_weight"
    }),
    <!-- Plot.dot(summaryMod4_tt_zR_OJS_T_selected, { -->
    <!--   x: "t", -->
    <!--   y: "mean_gr", -->
    <!--   fill: "green" -->
    <!-- }), -->
    Plot.dot(summaryMod4_tt_zR_OJS_T_selected, {
      x: "t",
      y: "weightDATA_std",
      fill: "orange"
    }),
    Plot.ruleX(summaryMod4_tt_zR_OJS_T_selected, {
      x: "first",
      y: 3,
      "stroke": "green"
    }),
    Plot.ruleX(summaryMod4_tt_zR_OJS_T_selected, {
      x: "last",
      y: 3,
      "stroke": "red"
    }),
    Plot.text(summaryMod4_tt_zR_OJS_T_selected,
      Plot.selectFirst({
        x: 1,
        y: 4,
        frameAnchor: "top-left",
        text: (d) => "Cohort = " + d.cohort
      })
    ),
    Plot.text(summaryMod4_tt_zR_OJS_T_selected,
      Plot.selectFirst({
        x: 1,
        y: 3.5,
        frameAnchor: "top-left",
        text: (d) => "Residual = " + d.meanResid
      })
    )
  ],
  facet: {
    data: summaryMod4_tt_zR_OJS_T_selected,
    x: "ind"
  }
})

18.5.7 Plot survival

Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Code

Plot.plot({
  width: width,
  height: 350,
  inset: 10,
  color: {
    scheme: "greys"
  },
  x: { label: "Age/season combination" },
  y: { label: "Probability of survival" },
  marks: [
    Plot.frame(),
    Plot.dot(summaryMod4_tt_zR_OJS_T_selected, {
      x: "t",
      y: "pSurv"
    }),
    Plot.dot(summaryMod4_tt_zR_OJS_T_selected, {
      x: "t",
      y: "enc8",
      fill: "orange"
    }),
    Plot.line(summaryMod4_tt_zR_OJS_T_selected, {
      x: "t",
      y: "pSurv"
    }),
    Plot.ruleX(summaryMod4_tt_zR_OJS_T_selected, {
      x: "first",
      y: 1,
      "stroke": "green"
    }),
    Plot.ruleX(summaryMod4_tt_zR_OJS_T_selected, {
      x: "last",
      y: 1,
      "stroke": "red"
    })
    ,
    Plot.text(summaryMod4_tt_zR_OJS_T_selected, 
      Plot.selectFirst({
        x: 11,
        y: 1,
        text: d => d.cohort
      })
    )
  ],
  facet: {
    data: summaryMod4_tt_zR_OJS_T_selected,
    x: "ind"
  }
})

18.5.8 Size-dependent survival

Code

  summaryMod4_tt_growth_phi <- MCMCsummary(object = out_NN_OBmod4$mcmc, params = c("phiInt", "grInt"), round = 3) |> 
    rownames_to_column(var = "var") |> 
    mutate(
      yoy = as.numeric(str_match(var, "([0-9]+), ([0-9]+)")[,2]),
      season = as.numeric(str_match(var, "([0-9]+), ([0-9]+)")[,3])
    ) |> 
    separate_wider_delim(var, "[", names = c("var1", "var2")) |> 
    dplyr::select(mean, var1, yoy, season) |> 
    pivot_wider(names_from = var1,  values_from = mean)

  summaryMod4_tt_growth_phiBetas = MCMCsummary(object = out_NN_OBmod4$mcmc, params = c("phiBeta"), round = 3) |> 
    rownames_to_column(var = "var") |> 
    mutate(
      beta = str_match(var, "([0-9]+), ([0-9]+), ([0-9]+)")[,2],
      yoy = as.numeric(str_match(var, "([0-9]+), ([0-9]+), ([0-9]+)")[,3]),
      season = as.numeric(str_match(var, "([0-9]+), ([0-9]+), ([0-9]+)")[,4])
    ) |> 
    separate_wider_delim(var, "[", names = c("var1", "var2")) |> 
    dplyr::select(var1, mean, beta, yoy, season) |> 
    pivot_wider(names_from = beta, names_prefix = "beta_", values_from = mean)
  
  summaryMod4_tt_growth_phiParams <-  
    bind_cols(summaryMod4_tt_growth_phi, summaryMod4_tt_growth_phiBetas |> select(-c(yoy, season))) 
    # mutate(
    #   phiInt = logit(betaphiIntTOut)
    # )
  
  weights = expand.grid(weight = seq(-2,2,0.25), yoy = 1:2, season = 1:4) |> 
    left_join(summaryMod4_tt_growth_phiParams) |> 
    filter(!(yoy == 1 & season == 2)) |> 
    mutate(pSurv = ilogit(phiInt + beta_1*weight + beta_2*weight^2))

  ggplot(weights, aes(weight, pSurv, color = factor(season))) +
    geom_point() +
    geom_line() +
    #theme_publication() +
    facet_wrap(~yoy)

18.6 Output model summary data for Xiaowei

Code

  #write.csv(summary_tt, './models/cmrNN_OB/tt/dataOut/xiaowei/summary_tt.csv')
  
  #write.csv(summary_tt_z, './models/cmrNN_OB/tt/dataOut/xiaowei/summary_tt_z.csv')